function out = runOptimization(setupFilter,params)
%% For logging
if setupFilter.MEXon == 1
    if setupFilter.optim > 0 && setupFilter.numCPUs > 2
        now             = clock;
        nameOfLogFile   = ['logFile_month',num2str(now(2)),'_day',num2str(now(3))...
            ,'_hour',num2str(now(4)),'_min',num2str(now(5)),'.txt'];
        %diary(nameOfLogFile)
    end
end
%% The estimation
paramsValues           = struc2values(params,setupFilter.selectParams);
exitflag               = 0;
if setupFilter.optim == 1
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    [paramsOptValues,fval,exitflag] = fminsearch(@objectFuncQML_ForOptim,paramsValues,options,setupFilter);    
elseif setupFilter.optim == 2
        % For Nelder-Mead with transformed space for the parameters
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    setupFilter.transformParamsOn = 1;
    paramsValues    = transformParams(paramsValues,setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);
    [paramsOptValues,fval,exitflag] = fminsearch(@objectFuncQML_ForOptim,paramsValues,options,setupFilter);
    paramsOptValues = unTransformParams(paramsOptValues,setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);    
    setupFilter     = rmfield(setupFilter,'transformParamsOn');   
elseif setupFilter.optim == 3
    % Matlab's gradient-based optimizer: Does not work well!
    %options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
    %    'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX,...
    %    'GradObj','on','Hessian','user-supplied','Algorithm','trust-region-reflective');
    %[paramsOptValues,fval,exitflag,output,lambda,grad,hessian]= fmincon(@objectFuncQML_ForOptim,paramsValues,...
    %    [],[],[],[],setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues,[],...
    %    options,setupFilter);
    % Own implementation of the Newton-optimizer
    setupFilter.transformParamsOn = 1;
    paramsOptValues = newtonOptim(paramsValues,setupFilter,...
        setupFilter.MaxIter,setupFilter.TolFun,setupFilter.TolX);   
    setupFilter = rmfield(setupFilter,'transformParamsOn');

elseif setupFilter.optim == 4
    % the CMAES
    InsigmaValues      = setupFilter.InsigmaValues;          %The std deviation in the initial search distributions
    sigma              = setupFilter.cmaesSigma;             %The step size
    opts.SigmaMax      = setupFilter.cmaesSigmaMax;          %The maximal value for sigma
    opts.LBounds       = setupFilter.lowerBoundsValues;      %Lower bound for params
    opts.UBounds       = setupFilter.upperBoundsValues;      %Upper bound for params
    opts.MaxIter       = setupFilter.MaxIter*100;             %The maximum number of iterations
    opts.MaxFunEvals   = setupFilter.MaxEvals*100;            %The maximum number of function evaluations
    opts.PopSize       = setupFilter.cmaesPopSize;           %The population size
    opts.VerboseModulo = 1;                                  %Display results after every 10'th iteration
    opts.TolFun        = setupFilter.TolFun;                 %Function tolerance
    opts.TolX          = setupFilter.TolX;                   %Tolerance in the parameters
    opts.Plotting      = 'off';                              %Dislpay plotting or not
    opts.Saving        = 'off';                              %Saving results
    opts.numCPUs       = setupFilter.numCPUs;
    %paramsOptValuesCMEAS = cmaes_dsgeDisplay_mp(@objectFuncQML_ForOptim,paramsValues,sigma,InsigmaValues,opts,setupFilter);    
    %paramsOptValues = paramsOptValuesCMEAS;
    [xmin,fmin,counteval,stopflag,outhist,bestever] = cmaes_dsgeDisplay_mp_grendel(@objectFuncQML_ForOptim,paramsValues,sigma,InsigmaValues,opts,setupFilter);    
    paramsOptValues = bestever.x;
elseif setupFilter.optim == 5
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    [paramsOptValues] = fmincon(@objectFuncQML_ForOptim,paramsValues,[],[],[],[],...
        setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues,[],options,setupFilter);
else
    paramsOptValues = paramsValues;
end
paramsOpt = values2struct(paramsOptValues,setupFilter.selectParams);
[Q,errorMes,model,setupFilter,outFilter] = objectFuncQML(paramsOptValues,setupFilter);
if errorMes == 1
    out = [];
    return
end

% We always want to have the model fit for surveys computed
modelPosterior = model;
modelPosterior.inclSurveys = 1;
if model.inclSurveys == 0
   % Solve for expected short rates
   y_select = find(strcmp(model.labely,'$r_t$'));
   vectorMom3 = zeros(1,model.ne);
   maxMatShortRate = max(setupFilter.matConShortRate);
   [rx,rxx,rss,rxxx,rssx,rsss] = CondMoments_Mom3_log(model.gx,model.gxx,model.gss,...
        model.gxxx,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,y_select,maxMatShortRate,setupFilter.appOrder);

   modelPosterior.matConShortRate = setupFilter.matConShortRate;
   modelPosterior.r0     = repmat(model.by0(1,1),maxMatShortRate,1);
   modelPosterior.rx     = rx;
   modelPosterior.Rxx    = 1/2*reshape(rxx,length(setupFilter.matConShortRate),model.nx^2);
   modelPosterior.Rxxx   = 1/6*reshape(rxxx,length(setupFilter.matConShortRate),model.nx^3);
   modelPosterior.rss    = rss;
   modelPosterior.rsss   = rsss;
   modelPosterior.rssx   = rssx;
end   
version        = 1+1;     %1 for evaluate at xhat and no account of state uncertainty
                          %2 for evaluate at xhat and accounting for uncertainty 
dimy           = setupFilter.numMacro+length(setupFilter.maturities)+length(setupFilter.matConShortRate);                         
outFilter.yhat = CDKF_y_posterior(setupFilter.data(:,1:dimy)',modelPosterior,setupFilter,outFilter.xhat,outFilter.Smat,version);
outFilter.what = outFilter.yhat - setupFilter.data(:,1:dimy)';

% Add a simulated sample path to the reported output if not already present
if setupFilter.numSim == setupFilter.burnIn
    setupFilterSim = setupFilter;
    setupFilterSim.numSim = size(setupFilterSim.shocks,2);
    [QSim,~,~,setupFilterSim,outFilterSim] = objectFuncQML(paramsOptValues,setupFilterSim);
    outFilter.outSim                = outFilterSim.outSim;
    outFilter.modelMomentsSim       = outFilterSim.modelMoments;
    outFilter.yieldsAllSim          = outFilterSim.yieldsAllSim;
    outFilter.yieldsAllSimWithNoise = outFilterSim.yieldsAllSimWithNoise;
    outFilter.ySimAsInData          = outFilterSim.ySimAsInData;
    outFilter.ySimAsInDataWithNoise = outFilterSim.ySimAsInDataWithNoise;
    outFilter.QSim = QSim;
end

% Adding term premia for nominal
TP_version   = 2;
model        = getTermPremia_loadings(model,setupFilter,TP_version);
termPrem     = getTermPremia(model,outFilter.xhat);
termPremSim  = getTermPremia(model,outFilter.outSim.xAll_cu);
T            = size(termPrem.TP,2);
TPdecom      = zeros(length(model.by0),T,model.ne);
modelTmp     = model;
for i=1:model.ne
    % Version 1
    %xhatDecom = 0*outFilter.xhat;
    %xhatDecom(1:model.nx1,:) = outFilter.xhat(1:model.nx1,:);   %endogenous states
    %xhatDecom(model.nx1+i,:) = outFilter.xhat(model.nx1+i,:);   
    %tmp = getTermPremia(modelTmp,xhatDecom);
    %TPdecom(:,:,i) = tmp.TP;
    
    %Version 2 (should be more accurate than version 1)
    tmpAll = getTermPremia(modelTmp,outFilter.xhat);
    xhatDecom = outFilter.xhat;
    xhatDecom(model.nx1+i,:) = 0*outFilter.xhat(model.nx1+i,:);   %turn-off shock i
    tmp = getTermPremia(modelTmp,xhatDecom);
    TPdecom(:,:,i) = tmpAll.TP-tmp.TP;
end
termPrem.TPdecom = TPdecom;

% Risk-adjusted CS regression loadings
% Convert into forwards
fwdPrem = nan(size(termPrem.TP));
for i = 1:size(termPrem.TP,1)-1
     fwdPrem(i,:) = (i+1)*termPrem.TP(i+1,:) - i*termPrem.TP(i,:);
end
% TPxhat and hence fwdPrem are scaled by 40000 to get risk premia in
% annualized basis points in getTermPremia. 
% Annualized yields in are not scaled. Hence, we need to undo the scaling of TP 
% to get basis points, but we still need to preserve the annualization i.e. the multiplication
% by 4.
numBoots = 5000;
CS_riskAdj = CampbellSchillerRegression_riskAdjust(setupFilter.yieldsAll,...
     termPrem.TP'/40000*4,fwdPrem'/40000*4,setupFilter.CSregress_m,numBoots);
%% Output and standard errors
out.exitflag        = exitflag;
out.paramsOpt       = paramsOpt;
out.model           = model;
out.setupFilter     = setupFilter;
out.outFilter       = outFilter;
out.termPrem        = termPrem;
out.termPremSim     = termPremSim;
out.CS_riskAdj      = CS_riskAdj;
epsValue = 1D-5;
% Compute standard errors if possible
if any(paramsOptValues-epsValue< setupFilter.lowerBoundsValues) ...
	|| any(paramsOptValues+epsValue > setupFilter.upperBoundsValues)
    out.SEcanBeComputed = 0;
    return;
else
    out.SEcanBeComputed = 1;
    setupFilter.setupProj.dispOn = 0;
    if setupFilter.seOn == 1
        for i=5:7
            setupFilter.epsValue = 10.^(-i);
            tmp = num2str(setupFilter.epsValue);
            if setupFilter.epsValue == 1D-2
               name = 'paramsSE_eps2';
            elseif setupFilter.epsValue == 1D-3
               name = 'paramsSE_eps3';
            elseif setupFilter.epsValue == 1D-4
               name = 'paramsSE_eps4';
            elseif setupFilter.epsValue == 1D-5
                name = 'paramsSE_eps4';
            else
               name = ['paramsSE_eps',tmp(5)];
            end
            disp([name,['  epsValue = ',num2str(setupFilter.epsValue)]])
            [out.(name)] = getQMLse(paramsOptValues,setupFilter,outFilter);
        end
    end
end
if setupFilter.MEXon == 1
    if setupFilter.optim > 0
        %diary off
    end
end

end